Plasma  Theory  &  Simulation  Group  ELECTRONICS  RESEARCH  LABORATORY 

DOE  Contract  No.  DE-ASO3-76F0O034-0E-ATO3-76ET53O64  University  of  California 

ONR  Contract  No.  N00014-77-C-0578  Berkeley  CA  94720 

SUMMARY  OF  PROGRESS  FOR  SECOND  QUARTER  -  1980 

This  is  a  summary  of  highlights  of  progress  made  by  our  group,  for  use  by  DOE  and  ONR  con¬ 
tractors.  The  details  appear  in  our  Quarterly  Progress  Report. 

LOWER  HYBRID  DRIFT  INSTABILITY  SIMULATIONS  IN  2D  USING  EI0HAR.  Our  study  of  the  lower 
hybrid  drift  instability  has  been  expanded  to  2d,  to  consider  nonlocal  effects  and  dif¬ 
fusion  phenomena.  Our  earlier  Id  work  modeled  electrons  using  a  linear  susceptibility; 
the  present  approach  uses  fully  nonlinear  particle  electrons.  Ions  are  still  assumed 
to  be  unmagnetized.  Initial ization  algorithms  are  developed,  and  some  test  runs  are 
described. 

MAGNETIZED  MULTI-RING  INSTABILITIES.  Instabilities  in  simulation  plasmas  arise  from 
the  loading  of  particles  onto  "spokes  and  wheels"  in  velocity  space,  a  technique  which 
is  commonly  used  to  reduce  noise.  Maxwellian  distributions  may  be  constructed  from 
(1)  N  equally  weighted  rings  and  (2)  N  uniformly  spaced  rings.  Threshold  values 
of  bip;/uc;  are  given  for  1,  for  Maxwellian  envelope. 

NONLINEAR  QUASINEUTRAL  SIMULATIONS  OF  FIELD-REVERSED  PLASMAS.  This  2d  code  is  now 
working  well,  and  has  been  applied  to  field  reversed  ion  layers.  Unstable  modes  with 
m*  2,3,4,  and  5  have  been  observed,  through  the  end  of  exponential  growth. 

FLUCTUATIONS  AND  LANDAU  DAMPING.  Thermal  fluctuations  in  a  simulation  plasma  and  the 
associated  resonance  broadening  and  Landau  damping  were  worked  out  and  compared  with 
simulations. 

MULTI-BEAMING  INSTABILITIES.  A  theory  for  the  multi-beam  instability  associated  with 
quiet  start  simulations  of  the  lower  hybrid  drift  was  made. 

ELECTROSTATIC  SHOCKS  IN  THE  AURORAL  MAGNETOSPHERE.  This  is  a  guest  contribution  from 
colleagues  at  the  Space  Sciences  Laboratory  (whom  we  helped  get  started)  describing 
particle  simulations  relevant  to  the  observations  of  electrostatic  shocks  or  double 
layers  in  the  auroral  magnetosphere  on  the  S3“3  satellite.  The  simulations  are  both 
Idlv  and  1d3v  (magnetized). 

LINEARIZED  3D  HYBRID  SIMULATIONS  OF  FIELD-REVERSED  SYSTEMS.  The  RINGHYBRID  code  has 
been  applied  to  exponential  rigid  rotor  equilibria,  using  a  scheme  whereby  the  coh¬ 
erent  part  of  the  perturbation  is  periodically  reconstructed  and  the  random  part  dis¬ 
carded,  so  that  stochastic-orbit  effects  are  suppressed. 


Our  group  currently  consists  of  five  Ph.D.  students,  programmer  (i),  engineering  aid  (i) , 
and  ourselves.  We  presented  two  papers  at  the  Sherwood  Fusion  Theory  Meeting  April  23_25 
at  Tucson,  AZ.  We  presented  four  papers  at  the  Ninth  Conference  on  Numerical  Simulation 
of  Plasmas  June  30-July  2  at  Northwestern  University,  Evanston,  IL.  One  journal  article 
appeared.  M.  J.  Gerver,  from  MIT,  was  a  visitor  for  a  week.  ONR  Contract  renewal  proposal 
was  mai led  Apri 1  24. 
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Station  I 

PLASMA  THEORY  AND  SIMULATION 

A.  LOWER-HYBRID  DRIFT  INSTABILITY  SIMULATIONS  IN  2d  USING  EZOHAR 
Yu-Jiqan  Chen  (Prof.  C.  K.  Birdsall) 

The  study  of  the  lower-hybrid  drift  instability  has  been  expand¬ 
ed  from  one  dimensional  simulation  using  a  hybrid  ESI  code  to  two  dimensional 
simulation  which  allows  us  to  study  the  nonlocal  effects  and  the  diffusion 
phenomena  using  the  EZOHAR  code. 

Similar  to  the  one  dimensional  simulation,  a  slab  shape  is  as¬ 
sumed  and  the  equilibria  are  functions  of  x  only.  The  drift  velocity  is  in 
the  y  direction  and  the  magnetic  field  is  in  the  z  direction.  The  simula¬ 
tions  are  carried  out  in  the  x-y  plane.  In  our  previous  one  dimensional  simu¬ 
lations,  ions  were  treated  as  fully  nonlinear  particles  and  electrons  only 
appeared  in  a  linear  etectron  susceptibility.  Now  both  ions  and  electrons 
are  fully  nonlinear  particles;  therefore,  all  the  electron  nonlinearities  can 
be  studied.  We  still  restrict  our  simulation  to  a  zero-beta  plasma. 

(1)  Equilibrium 

Since  the  mode  frequency  of  the  lower-hybrid  drift  Instability  is 
much  higher  than  the  ion  cyclotron  frequency,  the  ions  are  assumed  to  be  unmag¬ 
netized.  The  ion  pressure  force  is  balanced  by  the  ambipolar  electric  force 
In  the  equilibrium.  From  the  Vlasov  equation,  the  ion  distribution  function 
is  then  only  a  function  of  energy  H^  where 

Hj  ■  mjV^/2  +  e$(x)  . 

Let  the  Ion  distribution  function  fj(Hj)  be  exponential,  as 


(D 


fjfHj)  -  Cj  exp  l  - 


±mjV  +  e$(x) 


Rewrite  Eq.  (2)  as  a  product. 


fj(Hj)  -  n j (x)gj (v)  . 


Here 

2 

9j  (v)  -  C.  exp  (“in j v  /2T. ) 

is  the  Maxwellian  distribution.  The  ion  density  profile  is 


e  <Mx)  -»(0) 

T1 

n  is  the  ion  density  at  x«0. 
o 

Similarly,  the  electron  equi 1 Ibrium  distribution 
is  a  function  of  two  electron  invariants:  energy, 

H  ■  m  v^/2  -  e$(x) 
e  e 

and  the  guiding  center  position 


(x) 


"  exp 


X  ■  x  -  v  /« 
y  ce 


Let  the  electron  distribution  function  ^e(He,X)  be  given  as 


f  (H  ,X)  -  F (X)  exp 

B  6 


(2) 

(3) 

(4) 

(5) 

function 

(6) 

(7) 

(8) 


For  a  small  electron  gyroradius  and  a  slowly  varying  ambipolar  field,  i.e. 
a#  the  electric  potential  can  be  expanded  around  the  guiding  center 
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position  as 


$(X)  +  Ax  + 


where 


Ax  ■  x  -  X  »  vy/“ce  (10) 

is  the  electron  displacement  from  its  guiding  center.  Therefore,  Eq.  (9) 
becomes 


♦  (x)  -  <p (X)  -  v  E(X)/iu 

y  ce 


Substitute  Eq.  (11)  into  Eq.  (8)  and  rewrite  Eq.  (8)  as 


fe(He,X) 


F(X)  exp  1  + 


M  (X)  g  (v) 
©  © 


meVg(X) +e*(X)  me 


where  N  (X)  is  the  electron  guiding  center  density  profile  and 

g<[(»)  -  Ceexp^--i^*  [vy*y£(X)]2n  (13) 

is  the  electron  drifting  Maxwellian  distribution. 

From  the  Poisson  equation, 

■  4we[nj(x)  -  n#(x)]  ,  (14) 

and  Eq.  (5),  n^ (x)  and  n#(x)  can  be  determined  by  choosing  an  appropriate 
E(x). 
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An  example  of  an  equilibrium  is  shown  in  Fig.  1  when  the  E  field  has  the 


E(x)  ■  -4E  tanh*(x/L)  sech^(x/L)  .  (15) 

o 

(2)  Simulation 

The  subroutine  PTCLESET  of  the  code  EZOHAR  has  been  modified 
to  allow  particle  velocities  to  be  loaded  directly  from  an  input  file  named 
VEL.  If  the  parameter  IW  ■  1 ,  PTCLESET  will  get  the  particle  v^  values  in 
order  out  of  VEL  and  use  the  random  Maxwellian  loader  for  v  .  For  IW-2, 

y 

PTCLESET  reads  v  in  order  from  VEL  and  randomly  loads  v  in  a  Maxwellian 
y  x 

distribution.  When  IW-3,  both  vx  and  vy  values  are  obtained  from  VEL.  If 

IW* 4 (default) ,  the  Maxwellian  loader  is  used  for  both  v  and  v  . 

x  y 

In  our  LHD I  simulations,  all  four  options  were  tested  when  VEL 
was  set  for  the  quiet  start  Maxwellian  distribution.  Parameters  used  are: 
N^-64,  Nw«32,  Ne-Mj  -32768,  m,/m„-100,  ■  1 ,  (vc/vt.) _ -0.9  and 


pe'  ce  ’  E  ti 'max 


L/Aj»50.  Simulation  results  agree  very  well  with  the  lower-hybrid  drift 
instability  linear  local  theory.  The  nonlinearity  of  the  instability  has 
been  studied.  Saturation  levels  of  the  wave  were  in  fair  agreement  with 
that  predicted  by  ion  trapping.  Detailed  nonlinear  results  will  be  given 
in  the  next  Quarterly  Progress  Report. 


B.  MAGNETIZED  MULT  I -RING  INSTAB  I LITES 

Jin  Soo  Kim  (Prof.  C.  K.  Birdsall,  Dr.  B.  I.  Cohen) 


The  normalized  distribution  function  for  N  rings  Is  the  sum  of 


5-functions: 

Vvvi>  •  £.  S7-s<Vv«>  s(vl> 

S-l  JLS 


fh  n 

where  a8  is  the  weighting  of  the  s  ring  with  £  a^-1.  The  projects 
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of  f  on  the  v  axis  is 
o  x 

fo(vx}  "  /fo(Vvl)  dvydvn 


N 


~  —  ■  ■  0(v 

>1  “  V2  15 


V  ) 
X 


IS 


(2) 


where  6(x)  is  a  step  function,  equal  to  unity  for  positive  argument  and  zero 
otherwise.  This  f  (v  )  <s  sketched  in  Fig.  1  for  four  equal  density  rings 
(without  infinities). 

The  dispersion  relations  for  one  and  N  equal  density  ion  rings 

were  given  in  the  previous  QPR,  for  k||»0.  Adding  cold  background  electrons, 

2 

with  each  ring  density  reflected  in  its  gives 


tip-  N  a,  dJ^(u  )  Ifl 

°(“'kL)  ■  « ♦  #-  E  E  Pr 

e  S-1  s  is  s 


(3) 


where  u  *kv  /fl  and  J  is  the  Bessel  function  of  the  first  kind  of  order  2,. 
is  1J.S  s  i 


Adam  Drobot  in  a  private  communication  tells  us  to  look  for  another  term  in¬ 


side  the  sum,  from  relativistic  theory,  with  denominator  (u>  -  ,  which 


may  not  be  ignorable.  We  will  look  into  this. 


Typical  roots  of  D(w,k)  »0  are  shown  in  Fig.  2  for  four  unequal 


density  rings,  for  total  (u>p/nc)  .■  100 .  This  is  incomplete,  especially  at 


“LH  ^“pi  *  I©*  w^ere  some  other  behavior  might  be  evident.  Also  needed  is 


complex  u  for  one  ring  at  the  same  total  density,  for  comparison. 

Mike  Gerver  (private  communication)  points  out  that  we  should 


probably  always  find  instability  at  small  w/k^,  Inside  the  Innermost  ring. 


We  see  root  a.  in  Fig.  2  at  Rew»0  (so-called  zero  frequency  mode)  and  root  h_ 
at  Reu/k^  %0. I^v^;  if  we  went  to  larger  k,  we  might  see  roots  with  phase 


FIG.  2  Complex  u>  and  real  k  are  shown  for  four  rings  totaling  (01^/0)^,).  ■  100. 
Gyroradil  are  0.875»  1  •  75 *  2.63,  3.50  with  u>p*55.9,  35*5,  7-85,  and 
0.718,  respectively,  and  to  ^  —  1.0.  Ring  velocities,  relative  to  vt> 
are  0.4  (roughly  along  branches  c,e,f,g),  0.8,  1.2  (roughly  along  b,d) 
and  1.6. 


velocities  even  smaller  (than  h) .  The  slowest  ring  has  a^  -0.875  and  a. 

(at  v{)  is  2.19;  hence,  the  ring  velocities  are  radial  lines  of  slopes 
0.4,  0.8,  1.2,  and  1.6. 

We  are  extending  this  study  to  larger  values  of  the  parameters 
and  to  nonuniform  envelopes  of  the  rings,  especially  to  Maxwellian;  one  ob¬ 
ject  is  to  identify  the  number  of  rings  needed  for  stability,  a  sore  ques¬ 
tion  in  simulations.  Two  particular  loadings  of  rings  are  being  studied 
for  simulating  Maxwell ians:  (a)  N  equally  weighted  rings  and  (2)  N  uni¬ 
formly  spaced  rings  in  velocity  space,  following  the  multi -beam  studies  of 

■k 

Gitomer  and  Adam. 

By  choosing  suitable  v  (say,  3.5v  .  )  to  assign  values  of 

xnigA  cn 

vLS“|  v  ,  the  weighting  of  rings  is 

-v? c/2vth  /  N  ~v2  /2v2 

S  ■  \s  «  15  7  E  \s  «  15  th  .  <5> 

/  s-1 


For  the  equally  weighted  case,  we  know  that  a^  -a^ 


-  a^  -  -  and  then 


the  radii  of  rings  (v^)  in  the  velocity  space  are  chosen  by  the  following 
condl tion: 


,  r'“  •vi/7h  .  /  f  •vi/2vth  .  . 

*  Jo  •  ^V7o  •  ,iv^ 


The  a  and  f(v  )  are  shown  in  Fig.  3  for  both  spacings. 

S  X 


k 

S.  J.  Gitomer  and  J.  C.  Adam,  "Multibeam  instability  in  a  Maxwellian  Simu¬ 
lation  Plasma",  Phy.  FI.  pp.  719-722,  May  1976. 


The  dispersion  relations  for  the  two  models  are: 


D(u),k^) 


u>2  u>2 .  N  «  .  dJ2(u  )  J.Q. 

Vs  .--lar 

e  i  s”  1  is  is  i 


(5) 


Dtw.k^) 


u>-  N 


N  ui2  «• 


ps 


S=*l  |  lm-a 


1  <Uf(ra  ,)  M,  (6) 

(siixl)  d(suil)  u  -  fcft. 


in  the  presence  of  cold  background  electrons. 

Primary  interest  in  simulation  is  in  knowing  whether  a  Max¬ 
wellian  envelope  of  N  rings  is  stable  or  not,  in  order  to  choose  N  large 
enough  in  setting  up  initial  thermal  distribution.  Hence,  D»0  of  Eq.  (5), 
equally  weighted  rings,  was  solved  (by  Niels  Otani  using  Newton's  method) 
up  to  the  point  (w,k)  of  instability.  The  values  of  (w,k)  where  the 

2  2 

first  complex  roots  occur  are  given  in  Pig.  4  in  w-k  space  with  (N.^./fl. ) 

2 

as  parameters.  In  Fig.  5  the  (u>p./Q.)  sufficient  to  make  N  rings  unstable 
is  shown.  Note  that  use  of  Eq.  (4)  forces  a  v  ;  for  examples,  N»32, 

J.iuaX 


v  * 2.64v.  and  N»64,  v,  -2.89v  . 
imax  t  imax  t 
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c. 


FIELD-REVERSED  PLASMA  SIMULATIONS,  QU AS  I  NEUTRAL ,  IN  2d 
Douglas  Named  (Dr.  Alex  Friedman,  Prof.  C.  K.  Birdsatl) 


An  equilibrium  loader  has  been  developed  which  produces  rigid- 
rotor  field-reversed  equilibria.  This  loader  has  been  used  in  conjuction 
with  our  two-dimensional  quasineutrai  code,  AQUARIUS.  The  code  has  been 
applied  to  field-reversed  ion  layers.  Unstable  m-2,3,1*,  and  5  modes  have 
been  observed.  For  appropriate  parameter  ranges,  the  regions  of  instability 
have  been  found  to  correspond  to  those  predicted  by  Lovelace. *  A  descrip¬ 
tion  of  the  equilibrium  loader,  as  well  as  a  complete  discussion  of  ion- 
layer  stability  results,  will  appear  in  the  next  Quarterly  Progress  Report. 

D.  ORBIT  AVERAGING 

Vincent  Thomas  (Prof.  C.  K.  Birdsall,  Dr.  B.  I.  Cohen) 


Work  is  continuing  in  an  effort  to  Implement  orbit  averaging  in 
a  one  dimensional  electrostatic  code,  ESI.  At  the  present  time,  the  method 
appears  to  work  when  <  a>ce  and  k  •  B»0.  An  example  was  given  in  the  pre¬ 
vious  QPR  of  the  cold  lower  hybrid  dispersion  curve  (including  finite  kAx 
corrections) . 

For  the  case  u  >u>  (including  the  case  u  *0)  no  satisfactory 
pe  ce  9  ce 

results  have  been  obtained.  The  simulation  exhibits  a  violent  Instability  in 
the  first  few  time  steps  in  which  field  quantities  and  kinetic  energies  typi¬ 
cally  increase  by  a  factor  of  'vIOOO.  The  instability  is  saturated  by  trap¬ 
ping  in  the  potential  of  mode  1  and  is  characterized  by  X^^L,  the  length 

of  the  system.  This  behavior  is  insensitive  to  w  AT.  All  u>  AT  tried  from 

pe  pe 

the  stability  limit  of  the  leapfrog  algorithm  to  ui^AT  >  1000  have  been  unstable 

* 

R.  V.  Lovelace,  "Precession  and  Kink  Motion  of  Long  Astron  Layers",  Phys. 
of  Fluids  22,  No.  k,  pp  708-717  0979). 
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For  u  AT  *  10,  the  scheme  requires  taking  microsteps  for  the  Ions  as  well 

p6 

as  for  the  electrons  so  that  up;At|  never  violates  the  leapfrog  stability 
requirement  for  the  microsteps. 

A  different  method  for  biasing  the  equations  of  motion  was  pre¬ 
sented  by  8.  I.  Cohen  in  the  QPR  in  which  the  biasing  is  introduced  by  a 
backwards  derivative  of  the  velocity.  This  algorithm  has  been  tried  with 
orbit  averaging.  Similar  behavior  to  simple  orbit  averaging  was  obtained, 
except  when  the  biasing  was  very  large  (in  the  notation  of  the  last  QPR, 
e<1)  so  that  all  of  the  particles  became  frozen  at  equilibrium  positions. 

At  the  present  time,  we  have  developed  no  model  which  predicts 
the  stability  of  the  orbit  averaging  algorithm  as  functions  of  the  para¬ 
meters  of  the  problem,  particularly  as  u  ,  w  are  changed. 

C6  pC 

E.*  FLUCTUATIONS  ANO  LANDAU  DAMPING 

Douglas  Harned  and  Yu-Jiuan  Chen  (Prof.  C.  K.  Birdsall) 

(1)  Theoretical  Analysis  of  the  Effect  of  Thermal  Fluctuations 

The  spectral  density  of  a  warm  plasma  may  be  obtained  from  the 
fluctuation-dissipation  theorem: 


<E*  >  (k.iu) 
8ir 


co 


(k,u>) 


(1) 


Here  we  have  concerned  ourselves  with  longitudinal  fluctuations  with  8*0. 
The  expression  for  e(k,u)  is  obtained  from  a  Vlasov  model 


This  is  material  presented  at  our  regular  theory  and  simulation  seminar  on 
June  6,  1975,  and  considered  to  be  of  general  and  lasting  Interest,  hence 
included  here. 


where  the  summation  is  over  species.  Considering  a  two-component,  electron- 

ion  plasma  in  thermal  equilibrium  (i.e.,  T  »T.**A  -A  =A  ,  Maxwellian  velo 

e  i  e  I  o 

city  distribution),  e  becomes 

e(k,u)  -  1  -  “TT  [«(Z  )  +  W(Z.)]  (3) 

where  2  *iu/a)  kA_  and  W{Z  )  is  the  Fried-Conte  dispersion  function.  To  sub- 
s  ps  0  s 

stitute  this  into  the  fluctuation-dissipation  theorem  one  uses  the  relation- 


(e')2  +  (e")2 


where  eSe'+ie".  Upon  substitution  we  obtain 
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where  the  usual  analysis  has  been  used  to  obtain  the  imaginary  part  of  the 
dispersion  function,  i.e., 
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Ordinarily  the  assumption  is  made  that  kAn  «1  so  that  Z  becomes  large  in  the 


dispersion  function.  An  asymptotic  expansion  can  then  be  used  to  obtain  the 
spectral  density.  However,  here  we  will  be  interested  in  larger  values  of 
kAp  so  that  it  will  be  possible  to  observe  resonance  broadening  due  to  fluc¬ 
tuations  as  well  as  the  associated  Landau  damping.  In  such  a  case  (for 
kAp  *  .3)  the  expansion  technique  is  not  valid. 

The  dispersion  function  may  be  written  in  the  form: 

dy  exp  +  i  yf  Z  exp  ■  (6) 

A  program  SPECTRAL  was  written  to  evaluate  numerically  the  integral  in  Eq. 

(6)  and  then  to  compute  the  spectral  density  for  a  given  kA^  using  Eq.  (5)* 

For  large  values  of  Z  (Z  > 20) ,  the  large  argument  expansion  was  used  to  eva¬ 
luate  W(Z) .  This  program  was  used  to  provide  theoretical  curves  (for  kAQ  from 
0.1  to  1.0)  for  the  spectral  density  for  use  in  comparison  with  simulation 

values.  The  peak  corresponds  to  u  ,  and  the  width  is  a  measure  of  at,  , 

real  imag’ 

the  damping.  We  were  interested  in  the  values  near  kAp*0.1*  because  these 
values  seemed  to  be  most  amenable  to  simulation. 

2.  Simulation  to  Obtain  Spectrum 

The  simulations  were  performed  by  loading  a  random  Maxwellian 
velocity  distribution  of  8000  particles.  A^/Ax>10  was  used.  To  obtain  the 
spectral  density,  we  used  results  from  the  post-processor  ZED.  We  found 
that  to  obtain  reasonable  results  it  was  necessary  to  use  a  large  number  of 
plasma  periods  (we  used  ^3007^).  For  smaller  numbers  (e.g.,  20)  the  finite 
time  length  of  the  run  produced  an  inherent  (minimum)  line  width  that  was 
larger  than  the  resonance  broadening  due  to  Landau  damping  which  we  wanted 
to  study,  particularly  for  small  kA^.  Our  results  show  that  the  resonance 
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Z  exp 
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FIG.  1  Theoretical  spectral 
included  (near  u  -  0) . 


0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8 


FIG.  2 


Simulation  result  for  spectral  density,  for  m.  -*■«  (no  ion  response), 
for  kA  ■  .4,  with  about  300  plasma  periods. 
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broadening  seen  in  the  simulations  agreed  fairly  well  with  that  predicted 
theoretically.  Figures  1  and  2  show  the  theoretical  and  simulation  spectral 
density  plots  for  kXp*0.4.  The  theoretical  result  shown  includes  the  ef¬ 
fect  of  the  ions,  while  in  our  simulations,  the  ions  were  considered  to  be 
an  infinitely  massive  neutralizing  background.  This  is  not  a  very  important 
effect  since  ion  fluctuations  are  negligible  at  frequencies  as  high  as  the 
plasma  frequency.  However,  this  is  the  reason  for  the  low  frequency  peak  in 
Fig.  1. 

(3)  Simulation  of  Landau  Damping 

Having  obtained  the  fluctuation  spectra,  we  then  attempted  to 
see  Landau  damping  of  a  wave  excited  above  the  level  of  thermal  noise.  This 
requires  a  relatively  large  excitation  since,  for  the  number  of  particles 
ordinarily  used  in  simulations,  the  noise  level  is  high.  However,  with 
large  excitations  (at  levels  large  enough  to  let  particles  cross)  we  were 
able  to  see  evidence  of  Landau  damping.  The  damping  rate  for  kXD«  1  can  be 
obtained  by  measuring  the  width  of  the  Langmuir  resonance  and  by  direct 
measurement.  Figure  3  shows  the  result  for  kXp-0.4  and  Fig.  4  shows  the 
result  for  kXp»0.5.  The  agreement  between  simulation  and  theory  is  good 
as  shown  in  Table  1 . 

One  concern  was  that  for  such  large  excitations,  trapping  effects 
could  occur  which  might  obscure  the  damping.  In  order  for  trapping  to 
dominate  over  damping  one  should  have  >  |y|  where  y  is  the  damping  rate 
and  u>t  the  trapping  frequency,  ■  k>4$/m.  For  our  simulations  the  trapping 
frequency  and  the  damping  rate  were  comparable  at  the  initial  level  of  exci¬ 
tation.  However,  this  high  level  did  not  present  a  problem  because  after  one 
oscillation,  the  wave  had  damped  sufficiently  to  be  out  of  the  trapping  regime 
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Theoretical  and  simulation  values  for  Langmuir  frequencies 
and  damping  rates.  The  damping  rates  calculated  from  spec¬ 
tral  density  curves  were  obtained  by  measuring  the  width 
of  the  Langmuir  resonance. 


This  is  expected  since  any  large  amplitude  wave  that  is  constantly  damped 
has  particles  changing  from  trapped  to  untrapped  until  the  trapping  effect 
becomes  smallerthan  the  damping  effect.  In  Fig.  3  trapping  should  only  dom 
inate  during  the  first  oscillation.  This  may  account  for  the  different  ini 
tial  slope. 

(Langdon  comments  that  the  design  for  seeing  Landau  damping  is 


where 


kAv 


v  i 


<  kv  s  <  y,  , 

trap  trap  'Landau 

s  the  beam  spacing.) 

(The  work  reported  herein  is  presented  as 


interesting  but  not 


as  exhaustive  or  complete.) 
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F.  MULTI-SEAMING  INSTABILITIES 

Yu-Jiuan  Chen  (Prof.  C.  K.  Blrdsall) 

Multi-beam  instabilities,  MBI,  tend  to  plague  attempts  to  make 
highly  ordered,  hence  noiseless  (quiet)  initial  conditions.  Our  last  report 
on  these  problems  was  in  QPRI,  1979,  pp.  28,  29.  However,  they  have  con¬ 
tinued  to  plague  the  LHOI  work  and  have  been  observed  and  analyzed  in  our 
LHDI  studies  for  the  past  two  years;  see  especially  QPR  II,  1979,  Fig.  1  on 
p.  5  for  growth  rates  and  Fig.  2  for  LHOI  growth  plus  multi-beam  growth,  with 
analyses  in  the  text;  see  also  QPR  ill,  1979,  Fig.  2  on  p.  4  for  more  mixed 
mode  growth,  and  QPR  IV,  1979,  Figs.  2a, b  on  pp.  5,6  for  similar  growths. 

This  report  provides  theory  for  the  MBI  complementing  that  of 

•fe  ** 

Dawson  and  Gitomer  and  Adam,  and  was  given  as  part  of  our  plasma  theory 
and  simulation  seminar  June  5,  1979. 

Any  quiet  start  technique  allows  plasma  instabilities  of  interest 
to  be  followed  over  many  decades  of  growth  from  computer  roundoff  levels  to 
saturation.  However,  the  quiet  start  may  not  be  useful  for  weakly  unstable 
plasmas  or  for  wave  damping  studies  if  it  also  produces  the  multibeam  in¬ 


stab  i  1  i  ty. 


In  ESI  code,  there  are  two  ways  to  initialize  a  Maxwellian  velo¬ 


city  distribution,  random  or  quiet.  Here,  we  are  interested  in  the  quiet 
start  method  only.  In  this  method,  we  first  choose  the  number  of  beams 
NGR*N/NLG.  NGR  particles  are  uniformly  distributed  in  x«  (0,L/NLG)  where 
L  is  the  system  length.  A  set  of  discrete  velocities  are  found  from  the 


John  M.  Oawson,  "Plasma  Oscillations  of  a  Large  Number  of  Electron  Beams", 
Phys.  Rev.  1 1 8 ,  No.  2,  p.  381,  I960. 

S.  J.  Gitomer  and  J.  C.  Adam,  "Multibeam  Instability  in  a  Maxwellian  Simu¬ 
lation  Plasma",  Phys.  FI.  No.  5,  p.  719,  1976. 
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cumulative  probability  density  function 
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by  selecting  NGR  values  of  P  uniformly  in  the  interval  (0 . ,  1 . )  and  then 
finding  which  v(P)  corresponds.  For  a  Maxwellian,  ESI  uses 
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These  particles'  velocites  are  put  into  x  space  (L/NLG,L)  and  then  scrambled 
in  positions.  All  the  particles  have  equal  weights,  but  the  beams  set  up 
using  Eq.  (1)  are  unequally  spaced  in  velocity. 

Following  Dawson's  paper,  the  linearized  equations  of  motion  for 
these  beams  are 
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+  N . 
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where  n j  and  Vj  are  the  perturbations  in  the  density  and  velocity  of  the 
j1*1  beam  while  Nj  and  V.  are  the  corresponding  unperturbed  quantities. 
Assuming  that  solutions  have  the  form  exp[i  (kx-<ut)] ,  Eqs.  (5-7)  can  be  re¬ 
written  as 


i(u»-kVj)Vj  =  -eE/m 
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From  Eqs.  (8-10)  it  is  easy  to  obtain  the  well  known  dispersion  relation. 
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Equation  (11)  was  solved  using  SOLVER  ,  where  the  Vj  are  deter¬ 
mined  by  Eq.  (l).  For  NLG  ■  1  (meaning  one  particle  per  beam),  v  “0.141421, 
N-4096,  8192  and  16384,  the  real  part  of  the  frequency  of  the  most  unstable 
mode  (which  is  approximately  kV^^,  where  V|^  is  the  largest  velocity,  here 
taken  to  be  5vt),  and  the  corresponding  growth  rate  are  plotted  as  a  function 
of  kXp  in  Fig.  1.  Although  other  unstable  modes  are  not  shown  in  Fig.  1, 
all  of  them  have  similar  frequencies  and  growth  rates.  Figure  1  shows  that 

the  multi  beam  growth  rate  reduces  by  a  factor  of  N  at  large  k.  The  theore- 

H.  Stephen  Au-Yeung  and  Alex  Friedman,  "SOLVER:  An  Analytic  Function  Root 
Solving  and  Plotting  Package",  ERL  Memo  No.  UCB/ERL  M79/55,  Aug.  31,  1979- 
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tical  Langmuir  dispersion  curves  are  also  shown  in  Fig.  1.  As  0.3<kXp<0.6 
is  the  interesting  region  to  study  the  Landau  damping  of  Langmuir  waves.  Fig 
1  indicates  that  most  of  unstable  multibeam  frequencies  are  close  to  the 
Langmuir  frequency,  which  makes  it  very  difficult  to  excite  only  a  Langmuir 
wave  below  the  saturation  level  of  the  multibeam  instability.  A  similar 
dispersion  curve  for  N  *  16384,  NLG*128,  and  v^-0.049  was  given  by 
Gitomer  and  Adam. 

G.  ELECTROSTATIC  SHOCKS  IN  THE  AURORAL  MAGNETOSPHERE 

Mary  K.  Hudson  and  Oouglas  W.  Potter  (Space  Sciences  Laboratory, 

UC  Berkeley,  Guest  Contributors) 

In  this  section  we  summarize  our  recent  plasma  particle  simula¬ 
tions  relevant  to  the  observations  of  electrostatic  shocks  or  double  layers 
in  the  auroral  magnetosphere  on  the  S3~3  satellite.  First  we  review  our  one 
dimensional  in  space  and  velocity  (idlV),  unmagnetized  simulations  for  which 
we  used  an  electrostatic  particle  code  ESI  (Birdsall  and  Langdon,  1978). 
These  simulations  were  useful  for  learning  how  to  use  the  code  and  for  check 
ing  that  the  results  were  consistent  with  theory,  experiment ,  and  previous 
simulations.  We  then  modified  the  code  to  do  ld3V  simulations  in  a  magnet¬ 
ized  plasma  appropriate  for  the  magnetosphere.  Our  initial  results  are  re¬ 
ported  here  (Hudson  and  Potter,  1979). 

Table  I  summarizes  ours  and  other  double  layer  simulations  to 
date.  Most  previous  work  has  been  done  with  no  magnetic  field  imposed.  Our 
1d3V  simulations  consist  of  infinite  plane  charge  sheets  at  an  angle  ir/2  +  0 
with  respect  to  the  magnetic  field  which  is  in  the  Z  direction,  as  shown  in 
Fig.  1;  the  simulation  axis  of  charge  sheet  motion  is  along  X.  The  charge 
sheets  flow  along  B  when  an  electric  field  or  current  is  imposed  in  the  Z 


Solutions  for  the  most  unstable  mode  multibeaming  instability  with 
N- 4096(A),  8192 (x) ,  and  l6384(#)  beams.  Sketched  in  are  the  Landau 
roots  labelled  L. 
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Table  I 
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Goertz  &  Joyce 
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Hubbard  &  Joyce 
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Eo 
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de  Groot  et  al . 

0,  B 
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1 ,2d3Y 

Sato  &  Okuda 

0 
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<Ae 

IdlV 

Hudson  &  Potter 

0,  B 

P 

m 

o 

o 

0 

ld3V 

direction,  and  the  charge  sheets  oscillate  about  B  at  the  electron  and  ion 
cyclotron  frequency  respectively.  Thus,  our  simulation  retains  the  dynamics 
of  current  flow  along  B  and  magnetized  electrons  and  ions,  while  reducing 
the  number  of  particles  to  be  followed  and  the  number  of  spatial  grid  points 
where  fields  are  calculated  cubicly  from  the  three  dimensional  case. 

Our  initial  runs  were  1  d  1 V  with  B*0  and  with  periodic  boundary 
conditions  imposed.  We  chose  an  ion-to-electron  mass  ratio  M/m* 16  and  ion 
temperature  T.  *0,  later  changing  to  M/m *100  and  T.  »T  ,  as  described  below. 
We  imposed  a  constant  electric  field  Eq  greater  than  the  runaway  electric 
electric  field  which  accelerates  electrons  to  their  thermal  velocity  in  one 
collision  time.  We  chose  Eq*.01  in  normalized  computer  units  summarized  in 
Table  II.  The  electron  charge-to-mass  ratio  was  one,  therefore  the  electron 
drift  velocity  V*E  t  should  have  increased  linearly  with  t  such  that  V  *  1  at 


t*100  (uipe*l).  The  initial  linear  increase  in  electron  drift  velocity  is 
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Table  II 

Normalized  Computer  Units 


4) 

Ae 

M/m 

Tl/Te 

t 

AX 

At_ 

IdlV 

1 

.25 

16 

0 

50 

81.9 

.4 

.2 

ld3V 

1 

.25 

100 

1 

128 

32 

.25 

.2 

L  3  system  length 
n  =  charge  sheet  density 
ax  =  space  grid  steps 
At  =  time  grid  steps 

However,  the  Buneman  instability  intervened  by  t *  100.  The  fre¬ 
quency  and  growth  rate  of  the  Buneman  instability  can  be  determined  from  the 
plot  of  fluctuating  E  /8ir  vs.  t  in  Fig.  3  (w“.025  and  y  m  .2  in  units  of 
oi  »1).  The  wavelength  determined  from  plots  of  density  fluctuations  vs.  X 
at  fixed  t  is  the  order  of  100Ap.  These  results  are  consistent  with  linear 
theory  (Buneman,  1959). 

The  linear  increase  in  electron  drift  with  time  is  apparent  in 
Fig.  2,  up  to  the  onset  time  of  the  Buneman  instability.  Current  saturation 
and  falloff  at  t ^  100  is  due  to  backscatteri ng  from  fluctuations.  This  is 
evident  in  plots  of  the  electron  distribution  function  at  fixed  times.  Fig¬ 
ure  4a  is  a  plot  of  the  electron  distribution  function  at  t»80.  The  freely 

accelerated  electron  Maxwellian  peaks  at  V--.8,  consistent  with  V-E  t,  and 

o 

there  are  few  electrons  near  V-0.  Figure  4b  shows  a  backscattered  electron 
population  at  t  ■  120  responsible  for  the  decrease  in  current  in  Fig.  2.  Eq 
continues  to  accelerate  the  electron  distribution  function  and  turbulence 
continues  to  decelerate  it.  A  freely  accelerated  electron  population  is 
apparent  at  V--2.  at  t  ■  200  and  at  V»-5.  at  t  ■  500  in  Figs.  4c  and  4d. 
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Figure  3 


Plot  of  22/8w  vs.  t  from  t  =  o  to  t  =  100  f*pe-l  ), 
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The  continued  increase  in  current  at  later  times  in  Fig.  2  is 
due  to  a  combination  of  current  carried  by  runaways  and  the  increase  in  the 
critical  drift  for  instability  V  * Ag,  the  electron  thermal  speed.  The 
latter  is  increasing  due  to  turbulent  heating  of  electrons  and  its  increase 
is  apparent  in  the  velocity  spread  of  the  electron  distribution  in  Fig. 

4a-d.  Thus  the  following  scenario  occurs:  when  V^>Vc>  instability  develops 
which  reduces  and  increases  ^c<x^e>  thereby  shutting  off  the  instability. 
However,  * EQt  increases  again  until  > V£  and  the  cycle  repeats.  The 
oscillatory  nature  of  the  instability  is  apparent  in  Fig.  5»  which  is  an 
extension  of  the  plot  of  electric  field  fluctuations  in  Fig.  3  to  later  times. 
This  oscillatory  behavior  has  been  reported  in  other  simulations,  cf.  David¬ 
son  et  al.  (1970)  and  the  review  by  Papadopoulos  (1977). 

In  our  1d3V  simulations,  we  have  imposed  an  electric  field  along 
B  with  the  simulation  axis  at  an  oblique  angle  6»ir/2-.1.  The  current  beha¬ 
vior  plotted  in  Fig.  6  is  similar  to  the  ldIV  case.  This  is  a  plot  of  * 
j  sin  .1  *  . 1jz>  and  is  therefore  reduced  by  a  factor  of  10  from  the  1 d 1 V 

case  for  the  same  E  .  Note  that  Buneman  saturation  of  current  occurs  slight- 
o 

ly  later  than  in  the  1  d  1 V  case  because  ions  are  now  warm  (T.  *T  and  M/m*  100 

i  e 

is  assumed)  and  contribute  Landau  damping,  which  raises  the  threshold  drift 

and  increases  time  t  ■  V/E  to  achieve  that  drift. 

o 

We  observe  oscillating  potential  structures  evolving  by  t  ■  100w  ^ 
at  Buneman  wavelengths  A  ■  85Xp  along  X.  This  is  the  dominant  mode  thus  far 
in  our  simulations  because  it  has  the  fastest  growth  rate  of  any  current 
driven  electrostatic  instability  for  currents  above  its  threshold  which  are 
produced  by  the  applied  Eq.  In  the  magnetospheric  case  of  weaker  currents 
other  modes  may  prevail.  We  plan  to  simulate  the  effect  of  other  current 


Figure  5 
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driven  modes  by  (1)  applying  a  weaker  electric  field  for  a  longer  time  to 
see  if  current  saturation  by  instabilities  with  a  lower  current  threshold, 
e.g.,  ion  cyclotron  or  ion  acoustic  in  the  presence  of  electron  heating, 
precludes  excitation  of  the  Buneman  instability  in  the  magnetized  case; 
and  (2)  by  applying  a  current  below  threshold  for  the  Buneman  instability. 

We  observe  the  coalescence  of  potential  fluctuations  into  a 
single  structure  by  t * 600w^  in  Fig.  7.  In  a  ldlV  unmagnetized  simula¬ 
tion  with  periodic  boundary  conditions  Biskamp  and  Chodura  (1973)  have 
shown  that  a  succession  of  wave  bursts  of  the  type  evident  in  Fig.  5  each 
doubles  the  wavelength  of  the  dominant  mode  by  coalescence  of  BGK  modes. 
The  dominant  mode  asymptotically  approaches  the  longest  wavelength  allowed 
by  the  size  of  the  simulation  system.  Figure  7  gives  the  impression  that 
if  periodic  boundary  conditions  were  not  imposed  which  force  the  potential 
to  be  the  same  on  both  sides,  a  single  potential  drop  would  exist  across 
the  system.  The  magnitude  of  e<J>/T  «*  2.4  x  10^  is  enough  to  accelerate  1eV 
electrons  to  keV  auroral  energies.  Note  that  this  electric  field  struc¬ 
ture  is  oblique  to  the  magnet ic  field  with  a  parallel  component,  analogous 
to  the  electrostatic  shocks  which  have  been  observed  on  the  S3~3  satellite 
(Mozer  et  al.,  1977).  We  plan  to  reduce  the  applied  current,  change  the 
boundary  conditions  from  periodic  to  nonperiodic,  and  determine  under 
what  circumstances  oblique  shocks  of  the  type  commonly  observed  on  S3-3 
occur  on  auroral  field  lines,  vs.  one  dimensional  double  layers,  where 
the  potential  drop  is  strictly  along  B. 


1  profile  vs.  X  at  t  -  600. 
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H.  FIELD  REVERSED  EQUILIBRIA,  INTRINSIC  STOCHASTIC ITY,  AND  COLLECTIVE 
STABILITY 

Alex  Friedman 


We  have  developed  a  computer  program,  RIG1DROTOR,  to  generate 
ax i symmetric  field  reversed  Vlasov  equilibria  for  use  with  our 
linearized  3d  hybrid  simulation  code  RINGHYBRID1 .  In  these  equilibria 
the  current  is  carried  entirely  by  hot  ions  with  gyroradii  which  may 
range  from  infinitesimal  (as  in  the  Hill  vortex  equilibria)  to  of  order 
the  system  size  (as  in  ion  ring  equilibria);  charge  neutrality  is 
assumed.  We  describe  the  equilibrium  calculation,  the  intrinsic 
stochastici ty  and  hence  local  instability  of  particle  orbits  in  these 
equilibria,  the  problem  this  instability  presents  for  studies  of 
collective  behavior,  and  a  means  of  partially  overcoming  the  difficulty. 

We  assume  an  exponential  rigid  rotor  equilibrium  distribution 
function  for  the  ions  of  the  form 

F  =  no(27nnT)-3/2  e“C(H  -  OP,)/T]t  (1; 


where  H  *  mv2/2  is  the  particle  kinetic  energy,  Pa  =  mrvs  +  qrA/c  is  the 
canonical  angular  momentum,  A  =  Aa(r,z)  is  the  azimuthal  component  of 
the  vector  potential,  T  is  the  temperature  in  energy  units,  and  n0  is 
the  density  on  axis.  The  velocity  distribution  is  Maxwellian  in  the 
rotating  frame,  and  the  mean  azimuthal  velocity  is  everywhere  rfl.  After 
integrating  vflF  over  velocities.  Ampere's  law  can  be  written: 

(VXVX A$)s  =  (47r/c)qrnn0eC<*rnA/c  *  ®r2n*/2]/T.  (2 


The  algorithm  used  to  solve  this  equation  is  similar  to  a  new  MHD 
equilibrium  algorithm2.  We  introduce  an  r-z  grid,  and  new  parameters  CT 
and  Ca  such  that 


T  *  CT 


“0  =  Cn  / 


where  *  the  value  of  the  flux  function  V  *  rA  at  the  o-point.  T 

and  n0  are  thus  known  only  after  the  solution  has  been  found.  We  make 


an  initial  guess  for 


compute  f 


T,  and 


compute  the 


corresponding  current;  invert  the  curlsquared  operator  for  an  improved 
A;  and  so  iterate,  applying  an  under-relaxation  factor  of  order  0.S  to 
h 

This  psper  was  presented  et  the  Ninth  Conference  on  Numerical  Simulation  of 
Plasmas,.  Narthw»«r«rr.  llnf,u^,  fyenatiy  m  -i»««- 
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the  current.  If  CT  and  Ca  are  not  introduced,  convergence  is  apparently 
not  achieved  for  all  solutions3. 

By  integrating  the  density  over  volume,  the  total  ion  charge  is 
obtained.  A  representation  of  the  distribution  function  is  obtained  by 
assigning  locations  and  velocities  to  particles  in  a  manner  which  yields 
the  correct  density  and  mean  azimuthal  velocity  in  each  grid  cell, 
within  errors  associated  with  finite  particle  number.  The  r-z  plane  is 
divided  into  small  subcells,  and  the  area-weighted  density  at  each 
subcell  center  is  used  to  compute  the  probability  that  the  subcell 
contains  a  particle.  An  isotropic  velocity  distribution  in  the  rotating 
frame  is  generated,  with  a  user-specified  velocity  cutoff  (for  any 
component)  of  up  to  six  times  v thermal* 

RIGIDROTOR  runs  on  the  NMFECC  CDC-7600  (and  is  available  through 
LIBRIS);  the  output  file  is  sent  to  the  CRAY-1  where  its  field  and 
particle  initial  condition  data  are  used  to  initialize  RINGHYBRID. 
RINGHYBRID  advances  the  particles  along  their  equilibrium  orbits  in  the 
time- independent ,  selfconsistent  zero  order  field.  For  stability 
studies  of  low  frequency,  nonaxi symmetr ic  modes  the  code  simultaneously 
advances  linearized  displacements  from  these  orbits,  as  well  as 
linearized  quasineutral  fluid  equations  for  cold  ion  and  electron 
components,  and  selfconsistent  first  order  E  and  B  fields,  which  are 
fully  3d,  When  first  order  E  and  B  are  set  to  zero,  the  linearized 
displacements  are  advanced  using  forces  arising  solely  from  the 
equi l ibrium  magnetic  field  and  its  gradients;  they  thus  serve  as 
diagnostics  of  the  local  stability  of  each  orbit. 

A  variety  of  field  reversed  mirror  and  ion  ring  equilibria  have 
been  obtained;  when  the  equilibrium  field  and  particle  initial  condition 
data  are  used  in  RINGHYBRID,  particle  moments  remain  nearly  constant  in 
time.  Orbits  are,  in  general,  quite  complicated,  and  many  (from  10%  to 
most)  are  stochastic4.  The  complexity  of  moderate-gyroradius  FRM  orbits 
arises  from  the  fact  that  particles  moving  along  closed  field  lines  may 
or  may  not  reflect  in  the  high  field  regions,  and  furthermore  may  pass 
near  the  field  null  during  their  gyratory  excursions.  Also,  the  rate  of 
a  particle's  azimuthal  drifting  motion  depends  stongly  upon  where  it 
lies  on  its  orbit  in  the  r-z  plane. 

When  an  orbit  is  stochastic,  the  separation  between  it  and  any 
neighboring  trajectory  increases  exponentially  in  time,  with  noise 
superposed.  In  the  nonlinear  2d3v  zero  order  behavior,  the  only 
manifestation  of  orbital  stochastic! ty  is  an  eventual  loss  of  left-right 
symmetry  due  to  the  exponential  divergence  of  "neighboring”  mirror  image 
trajectories  (due  to  roundoff  effects,  they  are  not  perfect  mirror 


< 
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images).  This  zero  order  stochastic  behavior  has  been  studied  both  with 
surface  of  section  plots  (r-r  at  midplane)  and  histories  of  lrl-r2|  for 
mirror  image  particles.  However,  the  linearized  simulation  calculates 
first  order  currents  from  the  separation  of  orbits  which  are  displaced 
from  each  other  by  an  infinitesimal  vector  for  all  time,  and  in  many 
cases  collective  behavior  can  be  masked  by  rapid  separation  of  orbits  in 
the  equilibrium  field.  The  stochastic  parts  of  the  displacements  are 
random  in  phase  (uncorrelated),  and  so  (if  there  were  an  infinitude  of 
particles)  they  should  make  no  net  contribution  to  any  macroscopic 
quantity  such  as  the  current  -  only  the  coherent  part  of  the 
displacement  associated  with  the  collective  mode  should  not  average  to 
zero.  With  a  relatively  small  number  of  particles  in  each  cell, 
cancellation  due  to  random  phases  does  not  occur,  and  so  in  practice  the 
current  at  any  given  time  tends  to  be  dominated  by  the  contribution  of 
the  one  or  several  most  unstable  trajectories.  Unless  the  collective 
mode  being  studied  has  a  higher  growthrate  than  the  fastest  orbit 
separation  rate,  the  collective  behavior  is  invisible,  and  the 
"single-particle  modes"  dominate  (other  particles  respond  to  the  fields 
arising  from  the  large  offending-particle  current). 

Attempts  to  reduce  the  problem  by  employing  larger  numbers  of 
particles  were  unsuccessful,  presumably  because  cancellation  of  growing 
exponentials  is  impossible;  also,  as  more  particles  were  loaded,  some 
tended  to  fall  on  trajectories  with  even  higher  associated  growthrates. 
In  a  few  marginal  cases  the  most  unstable  particles  can  be  removed  from 
the  computation  artificially,  but  often  too  many  particles  must  be 
excised.  Until  recently,  we  have  been  forced  to  choose  equilibria  with 
low  rates  of  orbit  separation,  and  this  has  proven  to  be  a  significant 
restriction  on  the  applicability  of  our  program.  We  have  found  one 
moderately  thick  ion  ring  equilibrium  which  does  not  exhibit  rapid 
stochastic  instability,  and  have  examined  precessional  and  kink  mode 
behavior  of  this  system  in  some  detail5. 

Most  recently,  we  have  developed  a  procedure  whereby  the  coherent 
part  of  the  first  order  distribution  is  periodically  reconstructed  and 
the  random  part  discarded;  the  procedure  requires  an  approximate 
knowledge  of  the  mode  structure.  As  a  specific  example,  we  consider  the 
tilting  mode  of  an  ion  ring;  this  mode  is  essentially  rigid  in  nature, 
the  mean  axial  displacement  of  particles  at  any  point  in  the  poloidal 
plane  being  Independent  of  position  in  that  plane.  Thus,  every  (say) 
half  gyroperiod  we  compute  the  mean  axial  displacement  and  first  order 
axial  velocity  of  all  particles  (complex  quantities,  since  the  tilt  may 
have  any  orientation  in  azimuth),  and  set  the  axial  displacement  and 
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velocity  of  each  particle  equal  to  these  mean  values;  radial  and 
azimuthal  quantities  are  set  to  zero.  This  effectively  enforces  a  rigid 
tilt  on  the  system,  and  sets  the  random  part  of  the  displacements  to 
zero.  The  reconstruction  must  be  done  often  enough  that  the  coherent 
part  always  dominates;  the  random  part  must  never  be  allowed  to  grow 
from  the  noise  to  a  large  value.  When  this  procedure  is  applied  we  are 
able  to  observe  collective  behavior  in  systems  where  previously  only  the 
stochastic  growth  was  evident.  Even  though  the  true  mode  may  not  be 
perfectly  rigid,  its  projection  onto  a  rigid  displacement  is  large 
enough  for  the  technique  to  work.  In  principle  one  might  alter  the  mode 
structure  variational ly  to  get  the  peak  (true)  growthrate. 

Future  reconstruction  schemes  might  average  globally  or  locally 
along  a  fieldline  (for  flute  or  ballooning  modes),  or  on  a  cell-by-cell 
basis;  the  latter  would  eliminate  the  requirement  that  we  know  the 
poloidal  mode  structure,  but  has  stiffer  statistical  requirements. 
Furthermore,  the  modes  are  still  restricted  to  be  fluidlike  in  that  the 
displacements  and  velocities  at  each  point  are  single  valued.  To 
properly  capture  all  resonant  particle  effects,  the  reconstruction  must 
be  done  in  the  full  phase  space,  and  some  knowledge  of  the  mode 
structure  is  required.  Other  possibilities  include  an  orbit  averaging 
scheme  (since  the  modes  in  question  are  generally  of  low  frequency, 
while  the  noise  is  often  of  higher  frequencies),  and  a  careful  spectral 
analysis  of  the  data  to  try  to  diagnose  a  small  coherent  signal  in  a 
noisy  background. 
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Section  II 

CODE  DEVELOPMENT  AND  MAINTENANCE 


A.  ESI  CODE 

See  "Orbit  Averaging"  in  Sec.  I,  Part  D.  See  "Mu I ti beaming 

Instabilities"  in  Sec.  I,  Part  F. 

B.  EMI  CODE 

No  special  progress. 

C.  EZOHAR  CODE 

See  "LHDI  Simulations  in  2d",  Sec.  I,  Part  A. 

0.  RINGHYBRID  C00E 

Alex  Friedman 

Applications  of  the  RINGHYBRID  code  to  exponential  rigid  rotor 
equilibria,  and  a  scheme  under  development  uhich  periodically  reconstructs 
the  coherent  part  of  the  perturbation  quantities  and  discards  the  random 
part,  are  described  in  the  abstract  of  a  paper  presented  at  the 
Ninth  Conference  on  Numerical  Simulation  of  Plasmas,  reproduced  elsewhere 
in  this  QPR.  Here  we  briefly  describe  other  developments. 

A  facility  for  use  of  the  TALLYC  postprocessor  has  been  added  to 
the  code.  Calls  to  the  BZQHAR  1 ibrary  routines  "timer*  and  "timend" 
cause  code  operation  to  be  interrupted  periodically,  the  program  counter 
value  being  noted  at  each  interrupt  ion.  The  data  are  written 
to  a  dump  file  colled  Kid run>0,  and  one  instructs  the  TALLYC  routine  to 
use  this  file  and  the  symbol  table  file  as  input.  Output  consists  of 
a  table  showing  how  much  time  is  used  in  each  subroutine,  and  (optionally) 
a  breaKdown  by  code  lines  of  the  time  used  within  selected  subroutines. 
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T8LLYC  and  SZ0H8R  are  available  from  rfilem  directory  . takeme  of  userno.  1234 
To  run  T8LLYC,  type  T8LLYC  l<idrun>0  x5050200  <subroutine  names>  /tv, 
uhere  x5050280  is  the  version  of  RINGHYBRID's  executable  file  used, 
containing  symbol  table  information;  T8LLYC  informs  the  user  of  the 
name  of  the  file  it  produces. 

The  user  can  refer  to  the  listing  of  RINGHYBRID  to  see  how  to  call 
"timer*  and  "t intend*  in  other  codes  -  *9’  is  the  logical  unit  number,  and 
"itally"  and  “2000b"  are  the  name  and  length  of  a  scratch  array  dimensioned 
in  the  cliche  "commons*.  "idrun"  and  "1*  are  the  header  and  number  of 
words  in  the  header  of  the  outputfile  to  be  produced  by  tally. 

(Jhen  this  facility  is  employed  in  RINGHYBRID,  we  find  that 
the  particle  mover  consumes  the  bulk  of  the  cpu  time,  but  that 
much  of  the  work  is  concentrated  in  the  gat her /scat ter  (indirect  address) 
parts  of  the  coding  which  perform  field  interpolat  ions  in  scalar  mode  - 
the  mover  itself  is  vectorized.  In  principle  it  is  possible  to 
voctorize  the  gat her /scat ter  also  in  this  code,  since  the  algorithm 
gathers  Re(Br) , Im(Br) , . . . , Im(Ez) ,8r (equil ibr ium) .BzCequ i 1 ibr ium)  •  14 
quantities  at  each  of  the  nine  grid  points  intersected  by  a  particle, 
and  scatters  Re(rho) . Im(rho) »Re( Jr) ,  Im( Jr) . . . . , Im(Jz)  ■  8 
quantities.  Thus,  minimum  vector  lengths  of  8  are  possible,  and  may 
be  implemented  at  some  point  in  the  future.  8  greater  improvement  might 
be  achieved  by  interleaving  the  arrays  so  that  all  rho.J  components 
at  one  cell  appeared,  then  all  rho.J  components  at  the  next  cell,  and 
so  on.  then  all  B.E.BCequil ibr ium)  components  at  the  first  cell.  etc.  8s 
particles  are  3  cells  wide  in  each  direction,  this  would  permit  a  minimum 
vector  length  of  24  in  the  scatter  operation  (and  42  in  the  gather). 
Unfortunately,  this  change  would  necessitate  a  total  rewrite  of  the 
code,  since  diagnostics  and  fieldsolving  are  based  around  the  present 
scheme,  wherein  Re(rho)  values  at  cells  l  through  (nz+2)*(nr+2)  ore 
stored  first,  then  Im(rho)  values.  Re(Jr)  values,  etc. 
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Also  added  was  a  facility  for  uriiing  a  file  mode<idrun>  uhich  cor  bo 
interpreted  by  the  postprocessor  ZED  (after  being  reformatted  by  a 
program  called  PASTEAF  and  sent  to  the  CDC-Z600) . 

Spectra  computed  by  ZED  based  on  components  of  the  mean  first  order 
displacements  (epsilons),  or  first  order  magnetic  field  in 
a  preselected  diagnostic  cell,  are  quite  clean.  This  work  is  in 
progress;  ue  have  enjoyed  a  col  1 aboral ion  with  Bill  Nevins  in  this  regard. 

A  side  benefit  of  the  cell -by-cell  reconstruct  ion  algorithm  being 
implemented  is  the  availablity  of  a  diagnostic  revealing  the  mode  structure 
(spatial  dependence  of  epsilon)  in  detail.  Lie  find  that  normal  modes 
such  as  the  ion  ring  tilt  are  not  perfectly  rigid,  but  rather  entail  displacement 
peaKs  at  various  points  in  the  poloidai  plane,  typically  at  largest 
radius  and  at  the  edges  of  the  ring.  Some  previous! y-run  problems 
of  interest  will  be  re-run  with  the  new  diagnostic. 

Finally,  an  extensively  revised  usage  manual  for  the  code  is  available 
from  filem  directory  .takeme  of  usernumber  1234  (on  the  7600).  New 
cade  variables,  and  revised  instructions  for  running  the  code,  are 
described.  Detailed  examples  have  been  added. 


Mm 
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E.  POLARES:  A  TWO-DIMENSIONAL  R-8  ELECTROSTATIC  SIMULATION  CODE 
Niels  F.  Otani  (Prof.  C.  K.  Blrdsall) 

We  have  vectorized  the  charge-weighting  routine  and  the  mover. 
The  code  now  runs  approximately  three  times  faster  than  before.  About  50 
usee  are  required  per  particle  per  timestep  for  24  azimuthal  Fourier  modes 
—  still  fairly  slow.  Vector ization  was  variously  performed  with  respect 
to  either  Fourier  mode  number  or  particle  index. 

A  full-dynamics  electron-mover  scheme  is  now  implemented.  This 
will  allow  more  flexibility  in  the  types  of  problems  we  can  examine  with 


the  code. 


Improvements  are  continuing  to  be  made  in  the  code  diagnostics. 


So  far,  it  is  possible  to  plot  any  two  of  r,  9,  vr#  and  v0  against  one 
another;  also  kinetic  energy,  electric  field  energy,  and  total  energy  his¬ 
tory  plotting  capabilities  have  been  installed.  Potential  contour  plots 
and  mode  energy  histories  are  contemplated  for  the  near  future. 

The  difficulty  mentioned  in  the  previous  QPR  concerning  the 
incorrect  dipole  oscillation  period  has  been  corrected.  The  period  ob¬ 
served  in  the  simulations  is  now  in  excellent  agreement  with  theory.  How¬ 
ever,  some  damping  or  beating  effect  of  the  oscillations  has  been  observed. 
Also  energy  conservation  is  only  good  to  25$  which  is  unacceptable.  These 
problems  are  currently  under  examination. 

F.  WAVES:  A  UNIFORM  TWO-FLUID  ELECTROMAGNETIC  DISPERSION  RELATION  SOLVER 
Niels  F.  Otani  (Prof.  C.  K.  Birdsall) 

The  dispersion  relation  solver  appearing  the  last  QPR  has  been 
completely  rewritten.  The  Muller's  method  solver  has  been  replaced  by  New¬ 
ton's  method,  since  all  the  roots  appear  to  be  real.  This  seems  reasonable, 
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since  the  only  free  energy  in  the  system  is  the  difference  in  the  species  tem¬ 
peratures;  relaxation  of  this  kind  of  free  energy  should  not  be  possible  in 
a  fluid  code . 

Operation  of  the  code  is  fairly  straightforward.  Initially,  a 
thorough  search  is  made  for  roots.  initial  guesses  are  made  at  intervals 
small  compared  to  the  smallest  character istic  frequency  of  the  system.  The 
fraction  these  intervals  are  of  the  smallest  characteristic  frequency  is  ad¬ 
justable.  An  attempt  is  made  to  converge  each  of  these  guesses  to  a  root. 
Having  found  these  roots,  two  options  are  offered.  The  user  may  request 
the  relative  values  of  the  wave  quantities  characteristic  of  each  root.  The 
possibility  of  degeneracy  is  allowed  for  in  this  option:  two  orthogonal 
eigenmodes  will  be  generated  for  a  degenerate  root.  Alternatively,  either  the 
wave  number  k,  or  the  angle  between  k  and  Bq,  0,  or  both,  may  be  allowed  to 
vary.  WAVES  will  then  follow  the  values  of  the  roots,  and  plot  the  resulting 
dispersion  relation. 


I 
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G.  REPORT  PREPARATION  USING  TRIX/RED  AND  REDPP 

Alex  Friedman 

A  system  for  numbering  the  equations  in  reports  formatted  using 
this  editing  system  has  been  developed  and  made  available  to  us  by 
Art  Ualstead  of  LLL  (current  affiliation  UCLA).  The  abstract  prepared 
by  Alex  Friedman  for  the  Numerical  Simulation  Conference  (appearing 
elseuhere  in  this  QPR)  was  prepared  using  this  system.  The  RED 
source  file  for  this  abstract,  and  a  COSMOS  deck  which  sends  test 
copies  to  the  local  Versa tec  printer,  are  available  from  filem  directory 
.takeme  of  usernumber  1234  -  they  are  called  “simconf"  and  “swaps im“ 
respectively  (refer  to  QPR’s  II,  III  of  1979  for  details  on  generating 
typeset  output,  and  for  r*her  notes). 

Note  that  equations  .re  defined  (using  macro  "eqt“)  before  they 

I 

are  employed  (using  macro  "eqn"  -  the  last  argument  is  the  number 
of  lines  to  be  3et  aside  for  the  equation).  Macros,  and  basic  page 
layout,  are  defined  at  the  top  of  the  source  file.  Note  also  the 
usefulness  of  the  "quotient"  facility  "qt "  to  define  frequently  occurring 
or  complicated  expressions. 

Finally,  it  is  possible  to  pref  ix  "chapter"  numbers  to  the  equation 
numbers  to  get.  for  example,  (5.13)  inslad  of  just  (13).  To  do  this, 
change  "ns(dum, <num>) "  to  "ns(dum,<5.num>) "  in  the  macro  definition 
for  eqn. 

One  drawback  of  this  system  is  that  l ine  numbers  as  they  appear 
when  the  source  is  listed  on  the  Versatec  do  not  agree  with  those 
evident  when  using  TRIX;  the  latter  are  correct,  and  so  the  user 
must  mentally  alter  the  listed  numbers  by  one  when  referring  to  a 
Versatec  listing.  (One  of  the  line  throws  is  not  recognized  by 
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lines  43-44,  and  so  for  line  numbers  in  the  bulK  of  the  text  after 
thi3  point  one  has  to  add  one  to  the  number  the  Versatec  prints  out.) 

Another  example  of  this  type  of  formatting  is  an  excerpt  from 
Ualstead's  thesis,  called  "art's. memo"  and  also  available  from  the 
filem  directory  .taKeme.  The  user  is  encouraged  to  format  and  print 
this  section  as  well,  since  it  provides  a  more  complete  illustration 
of  the  capabilities  of  the  system. 

H.  FREX 

H.  Stephen  Au-Yeung 

FREX  has  been  updated  on  June  9,  I960.  The  neu  version  has 
the  follouing  changes: 

(1)  The  output  file  name  now  begins  with  "f3"  instead  of  "f0“  to 
match  the  TV80LIB  neu  file  naming  convention. 

(2)  The  size  of  the  output  file  will  nou  be  expanded  if  necessary 
therefore  the  SIZE  command  is  not  required  for  large  files 
anymore . 
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Suction  711 

SUMMARY  OF  REPORTS,  TALKS,  PUBLICATIONS,  VISITORS 


Abstracts  for  two  papers  presented  at  the  Sherwood  Fusion 
Theory  Meeting,  April  23-25  were  in  the  previous  QPR. 

Four  papers  were  presented  at  the  Ninth  Conference  on  Numer¬ 
ical  Simulation  of  Plasmas  ?t  Northwestern  University,  Evanston  1L,  June 
30  -  July  2,  as  follows: 


Paper 

PB-6 

OB-3 

OB-5 

PC-6 


V.  Thomas  and  C.  K.  Birdsall,  "Plasma  Hybrid  Oscillations  as 
Affected  by  Aliasing"  {much  as  in  prior  QPR's) 

Alex  Friedman,  "Field  Reversed  Equilibria,  Intrinsic  Stochas- 
ticity,  and  Collective  Stability"  (See  Sec.  1,  Part  H  for  abstract) 

V.  Thomas  and  B.  I.  Cohen,  "Orbit  Averaging  in  an  Electrostatic 
Code"  (much  as  in  prior  QPR's) 

Jin-Soo  Kim  and  C.  K.  Birdsall,  "Magnetized  Multi-Ring  Instab¬ 
ilities"  (much  as  in  this  QPR) 

The  following  paper  was  published: 


C.  K.  Birdsall  and  Neil  Maron,  "Plasma  Self-Heating  and  Satur¬ 
ation  Due  to  Numerical  Instabilities",  J.  Comp.  Phys.  36_,  pp. 
1-19,  June  1980. 

We  were  visited  by  M.  J.  Gerver  for  the  week  of  May  12. 
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